% GAUTIER LE BIHAN - 2020
% Replication files for "Shocks vs Menu Costs: Patterns of Price Rigidity in an Estimated Multi-Sector Menu-Cost Model?" Review of Economics and Statistics
%
% This code produces Figure O

clear all
close all
clear matrix
clc


load ..\..\Estimation_param\MS_produits_1sect_mean_CPlus\actual_moments_k
weight=actual_moments_k(:,12);
sum(weight)

load ..\..\Simulations_VC\MS_produits_kur_rhofixe_kurcor\stat4
param=stat4;

weight=actual_moments_k(:,12);
sum(weight);
p0=param(:,1);
mu_c=param(:,2);
sig_a=param(:,3);
rho_a=param(:,4);

param0=[p0 mu_c sig_a rho_a];

cd 
load ..\..\Estimation_param\MS_produits_1sect_mean_CPlus\param_mean
param1= param_mean
p1=tanh(param1(:,1));
mu_c1=param1(:,2);
sig_a1=exp(param1(:,3));
param1=[p1 mu_c1 sig_a1];
size(param1(p1~=0),1)


load ..\..\Estimation_param\MS_produits_3sect_Cplus\param_3sectors
param2=param_3sectors
p2=tanh(param2(:,1));
mu_c2=param2(:,2);
sig_a2=exp(param2(:,3));
param2=[p2 mu_c2 sig_a2];
size(param2(p2~=0),1)


load ..\..\Estimation_param\MS_produits_5lamsectors_CPlus\param_5sectors
param3=param_5sectors
p3=tanh(param3(:,1));
mu_c3=param3(:,2);
sig_a3=exp(param3(:,3));
param3=[p3 mu_c3 sig_a3];
size(param3(p3~=0),1)
 

FFF=figure(1);
max_p0=0.2;
x=p0(p1~=0,:);
y=p1(p1~=0,:);
subplot(1,3,1);
%scatter(x, y,sz);
scatter(x,y,'*','r' )

hold on 
max_p0=0.2;%P0;
x=p0(p2~=0,:);
y=p2(p2~=0,:);
subplot(1,3,1);
%scatter(x, y,sz);
scatter(x,y,'^', 'b')
hold off
hold on 
y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')

%HLB
pbaspect([1 1 1])

axis([0 max_p0 0 max_p0]);
xlabel('\lambda baseline') % x-axis label
ylabel('\lambda alt.') % y-axis label

max_p0=max(mu_c1(mu_c1~=0 & mu_c<0.2));%P0
x=mu_c(mu_c1~=0 & mu_c<0.2,:);
y=mu_c1(mu_c1~=0 & mu_c<0.2,:);
subplot(1,3,2);
%scatter(x, y,sz);
scatter(x,y,'*','r' )
hold on 
max_p0=max(mu_c2(mu_c1~=0 & mu_c<0.2));%P0
x=mu_c(mu_c2~=0 & mu_c<0.2,:);
y=mu_c2(mu_c2~=0 & mu_c<0.2,:);
%scatter(x, y,sz);
scatter(x,y,'^', 'b')
hold off
hold on 
y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')

%HLB
pbaspect([1 1 1]);

axis([0 max_p0 0 max_p0]);
xlabel('\mu baseline') % x-axis label
ylabel('\mu alt.') % y-axis label


max_p0=0.12;
x=sig_a(sig_a1~=1,:);
y=sig_a1(sig_a1~=1,:);
subplot(1,3,3);
%scatter(x, y,sz);
scatter(x,y,'*','r')
hold on 
x=sig_a(sig_a2~=1,:);
y=sig_a2(sig_a2~=1,:);
scatter(x,y,'^', 'b')
hold off
hold on 
y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')

hold off
axis([0 max_p0 0 max_p0])
%lsline
xlabel('\sigma baseline') % x-axis label
ylabel('\sigma alt.') % y-axis label

%HLB
pbaspect([1 1 1])

 orient(FFF,'landscape')
   print('..\..\figures\secondround_est_HLB.pdf','-dpdf','-fillpage')
 orient(FFF,'landscape')
 
  print('..\..\figures\secondround_est.pdf','-dpdf','-fillpage')



FFFb=figure(2);
max_p0=0.2;
x=p0(p1~=0,:);
y=p1(p1~=0,:);

subplot(3,3,1);
%scatter(x, y,sz);
scatter(x,y,'*','r' )
y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')

axis([0 max_p0 0 max_p0]);
xlabel('\lambda baseline') % x-axis label
ylabel('\lambda alt.') % y-axis label
%pbaspect([1 1 1])

subplot(3,3,4);
x=p0(p2~=0,:);
y=p2(p2~=0,:);
max_p0=0.2;
scatter(x,y,'o', 'b')

y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')
axis([0 max_p0 0 max_p0]);
xlabel('\lambda baseline') % x-axis label
ylabel('\lambda alt.') % y-axis label

subplot(3,3,7);
x=p0(p3~=0,:);
y=p3(p3~=0,:);

scatter(x,y,'^', 'k')

y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')


axis([0 max_p0 0 max_p0]);
xlabel('\lambda baseline') % x-axis label
ylabel('\lambda alt.') % y-axis label

%HLB
%pbaspect([1 1 1])
max_p0=0.2;%P0
x=mu_c(mu_c1~=0 & mu_c<0.2,:);
y=mu_c1(mu_c1~=0 & mu_c<0.2,:);
subplot(3,3,2);
%scatter(x, y,sz);
scatter(x,y,'*','r' )
y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')
title('Aggregate Sector: Mean moments')
%HLB
%pbaspect([1 1 1])

axis([0 max_p0 0 max_p0])
xlabel('\mu baseline') % x-axis label
ylabel('\mu alt.') % y-axis label

subplot(3,3,5);
max_p0=0.2;%P0
x=mu_c(mu_c2~=0 & mu_c<0.2,:);
y=mu_c2(mu_c2~=0 & mu_c<0.2,:);
%scatter(x, y,sz);
scatter(x,y,'o', 'b')

y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')
title('Aggregate Sector: Three-sector economy')


subplot(3,3,8);
max_p0=0.2;%P0
x=mu_c(mu_c3~=0 & mu_c<0.2,:);
y=mu_c3(mu_c3~=0 & mu_c<0.2,:);
%scatter(x, y,sz);
scatter(x,y,'^', 'k')
y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')
title('Five-sector economy')

%HLB
%pbaspect([1 1 1])

axis([0 max_p0 0 max_p0])
xlabel('\mu baseline') % x-axis label
ylabel('\mu alt.') % y-axis label
%HLB
%pbaspect([1 1 1])

max_p0=0.12;
x=sig_a(sig_a1~=1,:);
y=sig_a1(sig_a1~=1,:);
subplot(3,3,3);
%scatter(x, y,sz);
scatter(x,y,'*','r')
y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')

hold off
axis([0 max_p0 0 max_p0])
%lsline
xlabel('\sigma baseline') % x-axis label
ylabel('\sigma alt.') % y-axis label
%HLB
%pbaspect([1 1 1])

subplot(3,3,6)
x=sig_a(sig_a2~=1,:);
y=sig_a2(sig_a2~=1,:);
scatter(x,y,'o', 'b')

y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')
axis([0 max_p0 0 max_p0])
%lsline
xlabel('\sigma baseline') % x-axis label
ylabel('\sigma alt.') % y-axis label
%HLB

subplot(3,3,9)
x=sig_a(sig_a3~=1,:);
y=sig_a3(sig_a3~=1,:);
scatter(x,y,'^', 'k')

y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')

hold off
axis([0 max_p0 0 max_p0])
%lsline
xlabel('\sigma baseline') % x-axis label
ylabel('\sigma alt.') % y-axis label

%HLB
%pbaspect([1 1 1])



 orient(FFFb,'landscape')
   print('..\..\figures\secondround_est_HLB.pdf','-dpdf','-fillpage')
 orient(FFFb,'landscape')
 
  print('..\..\figures\secondround_est.pdf','-dpdf','-fillpage')

  
  

  
  %%%%%%%%%%%% REF 1 ROUND 2 : check products outside X=Y for sigma param
    figure (6)
x=sig_a2(sig_a3~=1,:);
y=sig_a3(sig_a3~=1,:);
scatter(x,y,'^', 'k')

y1=[0 max_p0];
x1=y1;
line(x1,y1,'Color','black','LineStyle','--')

hold off
axis([0 max_p0 0 max_p0])
%lsline
xlabel('\sigma 3sect ') % x-axis label
ylabel('\sigma 5sect.') % y-axis label

  sigma0=exp(param0(:,3))
   sigma1=param1(:,3)
   sigma2=param2(:,3)
   sigma3=param3(:,3)
  
  
  mat_sigma=[actual_moments_k(:,1) sigma0 sigma1 sigma2 sigma3 sigma0-sigma1 sigma0-sigma2 sigma0-sigma3];
  
  figure (3)
  hist(mat_sigma(:,6))
  
  figure (4)
  hist(mat_sigma((mat_sigma(:,7)>-0.2),7))
  
  figure (5)
  hist(mat_sigma((mat_sigma(:,8)>-0.2), 8))
  
  p1_1 = prctile(mat_sigma(:,6),1)
  p5_1 = prctile(mat_sigma(:,6),5)
  p10_1 = prctile(mat_sigma(:,6),10)
  
 p5_3 = prctile(mat_sigma((mat_sigma(:,8)>-0.15),8),5)
  p10_3 = prctile(mat_sigma((mat_sigma(:,8)>-0.15),8),10)
    p25_3 = prctile(mat_sigma((mat_sigma(:,8)>-0.15),8),25)
    p50_3 = prctile(mat_sigma((mat_sigma(:,8)>-0.15),8),50)

  
  p5_2 = prctile(mat_sigma((mat_sigma(:,7)>-0.15),7),5)
  p10_2 = prctile(mat_sigma((mat_sigma(:,7)>-0.15),7),10)
   p25_2 = prctile(mat_sigma((mat_sigma(:,7)>-0.15),7),25)
    p50_2 = prctile(mat_sigma((mat_sigma(:,7)>-0.15),7),50)
 
  
    mat_sigmab=[actual_moments_k(:,1:6) sigma0 sigma1 sigma2 sigma3 sigma0-sigma2 sigma0-sigma3];
 
    
    mat_sigmab3=mat_sigmab(sigma0-sigma3<-0.023,:);
    mean(mat_sigmab3)
    
        mat_sigmab2=mat_sigmab(sigma0-sigma2<-0.023,:);
    mean(mat_sigmab2)
    
    
    
    
    bins=42;
max=0.1;
min=-0.2;
%bins=80;
binValues_f=(min:(max-min)/bins:max);
[h] = hist(mat_sigma(:,7), binValues_f); 
counts1 = h/size(mat_sigma(:,7),1) ;
  bar(binValues_f, counts1, 'r', 'barwidth', 1);

  
      bins=42;
max=0.1;
min=-0.2;
%bins=80;
binValues_f=(min:(max-min)/bins:max);
[h] = hist(mat_sigma(:,8), binValues_f); 
counts1 = h/size(mat_sigma(:,8),1) ;
  bar(binValues_f, counts1, 'r', 'barwidth', 1);

  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mu0=param0(:,2)
   mu1=param1(:,2)
   mu2=param2(:,2)
   mu3=param3(:,2)
  
  
  mat_mu=[actual_moments_k(:,1) mu0 mu1 mu2 mu3 mu0-mu1 mu0-mu2 mu0-mu3];
  mat_mu3=mat_mu(sigma0-sigma3<-0.023,:);
   
      mat_sigmabmu=mat_sigmab(mu0-mu3>0.02,:);
 
         mat_mub=mat_mu(mu0-mu3>0.05,:);
  
    